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Abstract 


The Nondominated Sorting Genetic Algorithm (NSGA) is adapted and used to obtain 
multiobjective Pareto optimal solutions for three grades of nylon 6 being produced in an 
industrial semibatch reactor. The total reaction time and the concentration of undesirable 
cyclic dimer in the product, are taken as two individual objectives for minimization, while 
simultaneously requiring the attainment of design values of the final monomer conversion 
and for the number average chain length. Substantial improvements in the operation of the 
nylon 6 reactor are indicated by this study. The technique used is very general in nature 
and can be used for multiobjective optimization of other reactors. Good mathematical 
models accounting for all the physico-chemical aspects operative in a reactor (and which 
have been preferably tested on industrial data) are a prerequisite for such optimization 


studies. 



Chapter 1 


Introduction 

Most real world decision making design problems require simultaneous optimization of 
multiple objectives depending upon the nature of the system. Multiple objective 
optimization problems are conceptually different from single objective optimization 
problems. Single objective optimization problems always find the global optimal solution 
(maximum or minimum depending upon the requirement of a particular system). In case of 
multiobjective optimization, a set of solutions is obtained which are not necessarily the 
global best solutions if any of the objectives is considered individually, but are relatively 
better feasible optimal solutions depending upon the requirements of the designer if all the 
objectives are considered simultaneously. This set of solutions is called the Pareto-optimal 
or nondominated solutions. In this study, multiobjective optimization of an industrial nylon 
6 reactor is being carried out, and Pareto-sets are generated. These help to channelize the 
thinking of a decision maker in identifying suitable operating or design conditions. 

A number of studies have been reported on the optimal temperature and initiator or 
monomer addition histories (or profiles) for free radical and step-growth polymerizations 
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in batch, semibatch, or plug-flow reactors. These have been reviewed by Farber.^ In most 
of the studies discussed in the review a single objective function is used which is a 
weighted average of a few individual objectives as, for example, (a) concentration of 
unreacted monomer in the product, (b) concentration of undesirable side products, (c) 
reaction time, (d) deviation of the number average chain length (pn) and/or polydispersity 
index (PDI) firom the desired values, etc. This type of scalar optimization approach suffers 
from the drawback that the results are highly sensitive to the values of the weightage 
factors used. In addition, there is a chance of loosing some optimal solutions.^ In vector 
or multiobjective optimization, the objective function, I , is a vector, comprising of a few 

individual objective functions, F . The optimal solutions obtained using this approach help 
a designer to make better decisions. Since most optimization problems in polymer reaction 
engineering deal with different interesting and often conflicting objectives, multiobjective 
function optimization offers tremendous prospects for use in the optimal design and 
operation of these reactors. Some studies'**^ have already been reported on the 
optimization of polymerization reactors using this approach, but none of these involves 
industrial reactors. In this work, we have optimized the operation of an industrial nylon 6 
reactor. This reactor has been simulated by our group earlier, and a satisfactory 
mathematical model is already available, which can be used confidently for optimization 
purposes. 

Most of the earlier studies on the optimization of polymerization reactors involving 
single as well as multiple objective functions, use Pontryagin’s minimum principle’"^'* or 
some pattern search techniques to arrive dt optimal solutions. These traditional techniques 
usually require a good initial guess of the optimal solution, i.e., of the control variable 
histories. It is well known that Pontryagin’s technique is particularly sensitive to the choice 



of the initial guess. In fact, for complex, real-life problems (e.g., for methyl methacr\iate 
polymerization, as well as for nylon 6 polymerization,’^'** the ‘window’ in which the initial 
guess of the control variable must lie, is extremely narrow, and one can easily get into 
problems of lack of convergence. One often has to generate solutions of similar or easier 
problems before one gets an idea of good initial guesses.These techniques are, thus, not 
suitable for on-line applications. A new and extremely powerful search technique based on 
the mechanics of natural genetics and natural selection, called genetic algorithm (GA) is 
becoming popular in the literature. This technique does not require initial guesses of the 
control variable. This is a very robust technique and can be made to converge to the global 
optimum even in the presence of several local optima. Details of GA as well as of its 
several adaptations and extensions, are available in several books. A computer code 
(SGA - Simple Genetic Algorithm) which uses the basic algorithm is available now,^’ and 
this can be used as such or as a part of any larger code using any adaptation of this 
algorithm. In this work, however, we have used the nondominated sorting genetic 
algorithm (NSGA^*^) which is quite similar to the simple genetic algorithm (SGA) except 
for some differences in the working principle of the selection operator. Details of this 
adaptation are discussed later. 

Pareto sets for the industrial nylon 6 reactor are obtained using this technique. It is 
found that significant improvements (reduction in reaction time, reduction in the 
concentration of undesired side product) are indicated over the current operation of the 
reactor. We understand that some changes in the operating variables have been made in 
industry along the directions suggested by our studies, and considerable improvements 
have, indeed, been achieved. 
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Chapter 2 


Formulation 

The industrial reactor^®'*^ studied herein is shown schematically in Fig. 1. It is a 
jacketed vessel with a low speed anchor or ribbon agitator used for mixing the highly 
viscous reaction mass. The reaction mass (s-caprolactam, Ci, waier. W, and other inert 
additives like Ti 02 , etc.) is heated by condensing vapor in a jacket in -which the fluid is at 
temperature, Tj , (a constant value independent of time). As the poKinsrization takes place 
the temperature, T, inside the reactor goes above 22(fC and vaporizaticm of SKraprolactam 
and water occurs, resulting in a build-up of the pressure, p, in the regnan above the liquid 
reaction mass. The pressure inside the reactor, p(t), is maintained suer, as to conform to a 
desired history by manipulation of a control valve. The latter allows rae vapor mixture of 
nitrogen (inert, used initially above the liquid), monomer and water -na pass through to a 
condenser at a prescribed rate, Vt (mol of mixture/hr). The presssure histories used 
currently (all current values are referred to as ‘reference’ conditions nenoted by ref) for 
three different industrial runs producing different grades of nylon 6 (’csaing different initial 
concentrations of water, [W]o), are shown in Fig. 2 (cur\'es marked Tsri '’). In this diagram, 
the pressure and the time, t, have been nondimensionalized using 

n=(p-Po)/(P:nax.ra-PD) > 
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Condensing 
Vapor at Tj 


Condensate 


Fig. 1 : 



Schematic diagram of the industrial semibatch nylon 6 reactor 
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Fig. 2 ; Dimensionless pressure histories used currently to produce three grades of nylon 6, 
using different values of [W]o.^^ Optimal histories corresponding to points, O; , in the 
Pareto optimal solutions also shown. Arrows indicate the values of tf / t^f for the three 

mo - 
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Details of the model used (see Appendix 1), operating conditions, and results of 
simulation, are available in Ref 12. Five model parameters were ‘tuned’ using only one set 
(one value of [W]o ) of industrial data, and it was found that this tuned model gave results 
which compared satisfactorily with industrial data for the other two grades of polymer 
without requiring any retuning of the model-parameters. This indicated that aU the 
physico-chemical phenomena playing an important role in the reactor, are being correctly 
modeled. The model^’^ is used without any change for multifunctional optimization. 

The kinetic scheme for nylon 6 polymerization is given in Table 1.^'^ This scheme 
incorporates three important reactions, viz., ring opening, polycondensation and 
polyaddition, as well as two side reactions involving the cyclic dimer, Cz- Because of the 
non-availability of precise information about the rates of reactions involving the higher 
cyclic oligomers, C3, C4,..., these have not been incorporated in Table 1. Since the cyclic 
dimer constitutes a major share of the total cyclic oligomers in the reaction mass, this 
simplification of the kinetic scheme is justified. The ‘state’ of the reactor can be well 
described by a set of fifteen state variables, x (=[xi, X 2 , . . . , Equations for these 


can easily be written using mass and energy balances, and by obtaining moments of the 
polymer chain-length distribution, and are given in Ref 12 as well as in Appendix 1. In 
general, the state variable equations can be written in the form: 


dx 

dt 


-= fi(x,u) ; i= 1, 2, ... ,15 


( 2 ) 


where x and u are the vectors of the state and the control variables. 



1. Ring Opening 


C.+W^ 




■> 


ki' 


n 

K, 




2. Polycondensation 


k, 

: 

s„ + ^ + W; n, m = U,... 

k;=-^ 


3. Polyaddition 


S„+C 


i<- 


'k, 


^n+1’^ ~ 1,2,... 


4. Ring Opening Of Cyclic Dimer 


k. 


c, + w^ 




.s. 


5. Polyaddition Of Cyclic Dimer 


S„ + C; ^ , n — l.,2, . . . 


k, =r 
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contd. 



k. = A" exp(-E7RT) + a: exp(-E:/RT)i;([sJj = kf +k:l;i[s J 

n=l 


K,=exp[(AS,-AH,/T)/R],i = l,2,..,^ 


i 

A? 

(kg/mol-h) 

. E? 

(cal/mol) 

a; 

(kg'/ 

mol^-h) 

nr 

(cal/mol) 

(calymol) 

AS. 

(cal/ 

mol-K) 

1 


1.9880 

xlO^ 

4.3075 

xlO^ 

1.8806 

xlO'^ 

H.91S0 

xlCr" 

-7.88460 

2 


2.3271 

xlO^ 

1.2114 

xlO'° 

2.0670 

xlO^ 

-5.9458 

xlO-" 

-0.94374 

3 

2.8558 

xlO^ 

2.2845 

xlO^ 

1.6377 

xlO*° 

2.0107 

xlO^ 

-4.0438 

X 1 Cr 

-6.94570 

4 

8.5778 

xlO" 

4.2000 

xlO^ 

2.3307 

xlO^^ 

3.7400 

xlO'’ 

-9.6COO 

xlCr' 

-14.5200 

5 

2.5701 

xlO^ 

2.1300 

xlO^ 

3.0110 

xlO^ 

2.0400 

xlO^ 

-3.1691 

xlO"" 

-0.58265 


Table 1: Kinetic Scheme For Nylon 6 Polymerization-^^’ and 
corresponding rate parameters. 


9 































The ODEs (ordinary differential equations) in Eq. 2 can be integrated using the 
D02EJF subroutine of the NAG library for any given u(t) and initial conditions. This 


subroutine uses Gear’s technique for integrating sets of stin ODEs^®. The presence of 
some discontinuities and the stiffness of the ODEs makes it diScult to use a constant error 
tolerance (TOL) in the computer code, D02EJF. Provision w-as made in the algorithm for 
self-adjustment of the error tolerance between 10"^ and If the set of ODEs cannot 
be integrated for a certain value of TOL, the simulation package decreases the value of 
this parameter by a factor of 10 and integrates the set of differential equations again from 
the previous value of t (for which converged results have been obtained). This is 
attempted between upper and lower limits of TOL. This simulation package is combined 
with the NSGA^^ optimization code for performing multifunctional optimization. 

Two objective functions (both to be minimized) are considered in this study. The first 
involves the nondimensional reaction time, tf/tfjef . The total reaction time, tf, is determined 
as the time at which the monomer conversion, conv, reaches the desired value, convfj-ef 
(the value obtained in industry presently) as well as the degree of polymerization, pnj , 
reaches the desired value, pnjcf • The monomer conversion in a semibatch reactor at any 
time, t, is defined as 


conv = 1 - 


f[c.] +C, 

F.[C1 


( 3 ) 


where F and [Ci] represent the total mass of the liquid reaction mixture and the 
concentration ( mol/kg ) of the monomer at any time, respectively, and subscript, 0, 
represents the initial values. Ci is the (cumulative) amount (mol) of monomer that has 
vaporized until time t. 

The first of the above requirements (on conv) ensures that the monomer recycling 


load remains the same as at present, while the second requirement (on [in ) ensures that the 



physical properties of the product polymer (determined by pn primarily) are the same as 
those of the polymer currently produced. One may or may not be able to satisfy both of 
these two requirements at t = tf simultaneously. To account for this, the integration of the 
model equations (Eq. 2 ) is terminated at the point where conv = convf.ref .This denotes the 
point tf . The value, ptuf , of pn at this point is incorporated in the first objective function, 

Ii , in the form of a penalty function with weightage factor, Wi 

tf +w,(l-p„f /p„,^)- ( 4 ) 

Minimization of Ii leads to an increase in the production capacity of the industrial plant 
(through tf/tfjef), while simultaneously giving preference to solutions satisfying the 
requirement, pnj = Pn,ref (alongwith convf = convf^f) . 

The second objective function, I2, involves the nondimensional concentration, 
[C2]f/[C2]f,ref , of the undesirable cyclic dimer in the polymer produced. These cyclics cause 
problems in polymer processing^'^^ and are removed by hot water extraction, which is an 
energy consuming process. Minimization of [C2]i/[C2]f.ref , thus, leads to improvement of 
the product quality, as well as reduces extraction costs downstream. The penalty on the 
violation of pn.f = Pn,ref is also included in I2 with weightage factor, W2 . Thus, 

l 2 =[C 2 ]f/[C,Jf,^+w,(l-p„., ( 5 ) 

The use of penalty functions involving pnj in both the objective functions, and the use of 
the condition, convf = convfjef, to stop the integration of the state variable equations, 
ensures that we obtain optimal solutions where tf/tf^f and [C2]f/[C2]£ref are minimized, while 
ensuring convf = convtref as well as pn^ = Pnjef- 

The two objectives (Eqs. 4 and 5 ) are conflicting in nature and, therefore, provide an 
excellent opportunity for carrying out multifunctional optimization. In the industrial 
reactor, the vapor release rate history, VT(t), and the jacket fluid temperature, Tj (a 
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constant value), are used as the control (or optimizing) variables. These two constitute our 
control vector, u . Thus, our optimization problem can be written in mathematical form, as 


Min I = [ I, ,T ]" 

Tj.VtO) - L 1 ’ - J 

where 


(a) 


= tf ^ ' tf.rcf + W, (1 - (b) 

If =[C 3 f /[C 3 ],,^ + W 3 ( 1 -p 1 ;> /p,,,^)^i=l,2,...,Np (c) 

convf=convfjrf (d) 

mass and energy balance equations (e) 

any additional physical limits on the control variables (f) (6) 


The solution of the multiobjective optimization problem described in Eq. 6 is obtained 
using the nondominated sorting genetic algorithm (NSGA)^^ adapted in this study so as to 
apply for control variables which are continuous functions of time, t. Details of this 
adapted NSGA are described below, and the associated flowchart is given in Table 2. The 
numbers in the discussion below refer to the numbered boxes in Table 2. 

4. At generation number, Ng=0, a population having Np chromosomes (members) is 
generated. Each chromosome in this population carries the information of one digitized 
control variable history [set of vapor release rate, VT(t), values] and one control variable 
value (a single value of the jacket fluid temperature, Tj). We discretize our first control 
variable history, VT(t), in terms of Ng, equispaced points in 0<t<tfo (^fo> ^ 
estimate of tf , is to be supplied). The second control variable, Tj , is appended in any 
chromosome as an additional [(Nga+ 1)*] point. Thus, each of the Np chromosomes (called 
strings) comprises of a sequence of (Np + 1) numbers (called substrings). Each of these 
substrings, in turn, comprises of a set of Nar binary numbers (0 or 1). Each chromosome. 
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8. Initialize front number : 


No 


9. front = front + 1 



15 










Table 2: Flowchart Describing the Adapted NSGA Used in This Study. 
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therefore, has +i) x binary digits, (see box No. 4 of Table 2). The Nehr 

individual binaries in each of the Np chromosomes, are generated using a random number 
generation subroutine.^'"°'^’ 

5. The complete binary string (sequence of Nchr binaries) of the i^^ chromosome, when 
decoded into real numbers, Up j. or Tj'^and interpolated (mapped) between the upper (u < 

u ) and lower (u > u""") bounds of the control variables u, gives a digitized Up ^-history 

(a set of Nga real values), [Upj,Up 3 ,...,Upl.j^], representing a VT(t) history as well as a 

value, Tj'\ representing a single value of the jacket fluid temperature Tj, corresponding to 
that chromosome. Thus, there is a set of Np chromosomes, each representing a digitized 
up(t)[= V.i-(t)] history and a value of T;, each appropriately coded, in the form of a string 
of Nchr binaries. The minimum difference between the digitized values of the control 
variable (at any value, tu ) of two different chromosomes is (u”" - Uplf) / (2^" - 1). This 

is the accuracy to which any particular u? can be determined for different members of the 
population in a particular generation.. 

6. The decoded and adaptively mapped, discretized values, Upj., are curve-fitted 

piece- wise (splines) to obtain a continuous function, Up^(t). A piece-wise cubic Hermite 
subroutine (EOIBFF from the NAG library) is used to do this. This continuous function is 
again digitized to give Nsim(2:Nga) values of the control variable, [Up]; 1=1, 2,...,Nsim ]. 

7. These more closely spaced, discretized values of U^p ^ (t) and the value of Tj'^ are fed 
to the simulation package, D02EJF (of NAG library) which integrates the state variable 
equations (Eq.; 2 or Appendix 1) starting with the given initial conditions'^ and terminating 
at the stopping condition, convf = conv^f. The simulation program stores the values of 
each of the state variables, x'‘’(j), at every Hin, intermediate values of t, such that there are 



sets of X values until convf - conv^jef . These detailed histories could be printed out later 

for the optimal solutions. The values of the two objective functions I', ' and I. [at the final 
reaction time t = tf , see Eqs. 6 (b) and (c) ] are computed. One additional point needs to 
be emphasized. The computer codes involving GA usually maximize a fitness function, 
rather than minimize objective functions, m = 1, 2. Hence, we define fitness 
functions to convert the minimization problem to an equivalent maxiiuzation problem as 
follows; 


F<'>.l/(l+lS'>) 


(7) 


It may be added that the technique described till this stage is quite similar to that^ 
developed for the optimization of MMA polymerization. All the feasible points are 
identified (for plotting) at this stage. The feasible points or chromosomes are those for 
which conVf = convf,ref + 0.009 and (1- Un/tinjrf) = ± 0.009. The satisfi"action of the latter 
requirement means that the penalty value in Eq. 4 and 5 are negligible. 

9. A chromosome, ii , is said to be dominated by another chromosome. i 2 , (for the 
present problem of minimization of I or maximization of F), if 

(a) 

as well as 

Fj''^<Fj'-' (b:» 

then 

ii is dominated by i 2 (c) (8) 

For example, in Fig. 3 point 3 is dominated by point 5. We test each of the Np 
chromosomes in the population against all others to sort out all dominated ctiromosomes. 
As soon as a chromosome is found to be dominated, if is not checked for dominance with 



any other chromosome in ti:e population. When all chromosomes have been checked for 
dominance, and all domir.ated chromosomes have been identified, the rest of the 
chromosomes are given a front number, FRONT=l. These chromosomes having 
FRONT=l, are called nondominated chromosomes. In Fig. 3, points 4, 5, 6 will constitute 
the chromosomes assigned FRONT=l. All nondominated chromosomes are then assigned 
a dummy fitness value, F* equal to Np. Thereafter these dummy fitness values are modified 
according to the sharing procedure described in item 10 below, to assign a shared fitness 
value. Sharing is done to maintain diversity in the nondominated chromosomes.^^ In order 
to identify chromosomes for other fronts, we temporarily discard all nondominated 
chromosomes. The remaining chromosomes are again checked for dominance using Eq. 
8 and new nondominated chromosomes are sorted out and given a front number, 
FRONT=2. Again, the new nondominated chromosomes (in FRONT 2) are given a 
dummy fitness value, F,*, which is slightly smaller than the lowest of the shared fitness 
values of the previous front. The sharing of the dummy fitness values is performed again, 
and a shared fitness value is assigned to each nondominated chromosome. This procedure 
is continued until all Np chromosomes have been given a front number. 

10. Sharing : Sharing is performed among the members of the i*^ front (having n; 
members) using the following procedure: 

(a) For each chromosome, j, in front i, the dimensionless distance, djk , of this chromosome 


from any other chromosome, k , (including j) in the (same) front is calculated using 


d^, = z [(ug - ug)/! u^' - u”r)]' +[(x" -T^‘=y(Tr -ir)]' 


(b)We calculate the niche count , mj , using 
m. = ISh(d,) 


k=: 


( 10 ) 




Fig. 3 ; An example describing the concept of dominance . 
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where 



1 - 


" d, Y 




snare ' 

0, Otherwise. 




( 11 ) 


ashare is giveii by 


1 

*^share ~ 2q^’''<^'ga''W ' (^^) 

In Eq. 12, q is the number of Pareto-optimal solutions desired (we have used q = 15 in our 
study). The parameter, a, is an exponent which controls the sharing effect (we have used 
a = 2 in the present study). 

(c) The dummy fitness, , of each chromosome, j , in front, i , is modified by dividing F* 

/ 

by the chromosome’s niche count, mj , to calculate the shared fitness value, F,* , as follows; 


m. 


( 13 ) 


1 1 . The stochastic remainder roulette wheel selection is used on the shared fitness 
values, and a mating pool of N? chromosomes is generated. This procedure involves 
proportionate selection, where first, the number of copies made of each chromosome is 

t t t 

equal to the integer part of the value of F* / F,‘ . Here, F* is the average of the shared 

fitness values of all the N? chromosomes in the population. Additional copies of the j'*’' 
chromosome in the i*‘ front (to make a total of N? in the mating pool) are made 
thereafter, using a roulette wheel with probability proportional to the fractional part of 

tf 

F* / F* . 

u u . ' 

12. After the mating pool is created, crossover and mutation take place to produce the 
new population (next generation). These operations take place at the chromosome (binary) 



level. Two chromosomes are randomly selected from the mating pooh ^ crossing she is 
selected (randomly again), and portions of the chromosomes before and after the crossing 
site are exchanged. For example, for seven-bit chromosomes with crossing she after the 
third binary, the crossover is described by the following ; 

lOOlllll 100 0100 

► 

lioloioo 110 1111 (14) 

(old generation) ■ (new generation) 

While performning crossovers, only pcNp chromosomes are crossed, the remaining being 
left untouched (pc is referred to as the crossover probability). 

Another operation, called mutation, is also used to improve the next generation. The 
mutation operator changes a binary number from 1 to 0 or vice versa, with a probability , 
Pm. This operation is carried out for each of the NpN^j^ bits in the population, again using 
appropriate random numbers. The need for mutation leads to a local search around the 
current solution and helps maintain the diversity of the population. 

This completes one generation of NSGA. These sets of operations are carried out 
from one generation to the next generation until the number of generations equals the 
maximum number of generations specified at the starting of the program as an input 
parameter. 

In a nutshell, the feasible nondominated solutions are emphasized in every generation 
by the action of the reproduction operator of NSGA. Coexistence of multiple feasible 
nondominated solutions is ensured by using the sharing technique. Thereafter, crossover 
and mutation create new and, hopefully, better feasible nondominated solutions. This 
process continues and finally the best nondominated solutions of a population converge to 
the true Pareto-optimal solutions. Since diversity in the best feasible nondo m i n ated 



solutions is maintained in each generation, the final population is expected to capture a 
number of different Pareto-optimal solutions simultaneously in one tun. 



Chapter 3 


Results and Discussions 

Several checks have been made to ensure that the computer code was free of errors. 
For example, the values of conv, [C 2 ], Hn , etc. , for any one chromosome were found to 
be the same as those obtained using our earlier'^ simulation code [with the same VT(t) and 
Tj]. Also, Pareto optimal solutions for some simpler problems^ were obtained using our 
code, and they were found to be same as solutions obtained earlier. The CPU time ([W]o = 
3.45%) for 20 generations was found to be 144 min. on a mainframe HP8000S/950 
supermini computer. 

The upper and lower limits of each of the Nga control variables (for [W]o = 3.45%), 
Up“ and Up” , k = 1, 2 ,..., Ng, , as well as of the jacket fluid temperature, Tj““ and Tj"™ , 
were first determined. It was found that these (slightly restricted) search ‘windows’ for 
VTA'’T,ma.'c.ref and Tj were necessary, and that very large search domains on these two 
control variables led to severe computational problems with no Pareto solutions being 
obtained. Since our earlier multiobjective optimization study^* of this reactor (v/ith the 
shape of the pressure history fixed and described by four values only) gave us an optimal 
VtA/’w^f history and an optimal value of Tj , we selected the windows of these control 
variables around these previous solutions as 



(15) 


0.4 X ( Vt /\ T,n,ax.rcf ) ^ '^T /'^T.ma.^ta ^ 1-6 X ( V-j. / Vj j 

at corresponding values of t / tf . The window for Tj was taken as 220'^ C < X < 270’’ C. It 

j 

was found that the window for VTA'TJna.^u^f at t / tf = 0 had to be widened beyond that 
given in Eq. 15, before solutions could be obtained. The exact window used for 
VTA^T.ma.\.ref is givcu in Table 3, while Fig. 4 gives the window as well as the optimal values 
obtained earlier.^* Any further increase in the width of the window for VT/VTjniLx.ref reduces 
the number of Pareto points obtained (without any change in the optimal solutions), and 
no further widening of the width of the window was attempted. The values of the 
computational parameters used in the optimization study are given in Table 4. The initial 
and other values associated with the reactor are the same as those given in Ref 12. 

Figs. 5-9 show the feasible solutions from generation to generation for [W]o = 3.45%. 
Each of these points satisfy all the constraints of the system (Eqs. 6(d) - (f), and Pn = Pajef 
within error bounds). Convergence to the optimal Pareto set is observed to take about 8 
generations. The number of Pareto points obtained finally (ten) is found to be of the same 
order as q (15), as expected. The solution in Fig. 9 is referred to as the ‘reference’ run. 
The Pareto sets obtained for the other two water concentrations, i.e., [W]o = 2.52% and 
4.43%, are given in Figs. 10 and 1 1 . It may be added that the values of [C 2 ]fjef are different 
for the three cases, while tcnjf is the same. The values of Up“,Up“ (k = 1, 2, . . . , Nga), 
T™“ and T™” , as well as all the other computational parameters, for all these cases are 

taken to be the same as given in Table 3 and 4. 

After the generations of the Pareto sets, a designer can choose an operating point 
(preferred solution) from these sets using his intuition (or other information not included in 
Eq. 6). One technique which can be used is to select the preferred solution as a point 
nearest to a point, Ui , called utopia. The points, Ui , are obtained by performing two 



Fig. 4 : Dimensionless vapor release rate histories for [W]o = 3.45 %, for three points (tf /. 
tfjrf = 0.45, 0.54, 0.65) on the Pareto optimal set. and denote the lower and upper 
limits of the window used in this study, while denotes the optimal values of 


obtained in Ref. 28. 




k 







1 

0.05 

0.3 

0.01 

2 

0.06 

0.1 

0.02 

3 

0.2 

0.32 

0.08 

4 

0.3 

0.48 

0.12 

5 

0.2 

0.32 

0.08 

6 

0.07 

0.112 

0.028 

7 

0.03 

0.048 

0.012 

8 

0.05 

0.08 

0.02 

9 

0.02 

0.032 

0.008 

10 

0.02 

0.032 

0.008 

. 


* k •— 1,2,..., Nga 


Table 3: Window for Yj/ VT.n^^f Used for our Study [W]o = 3.45 

































( \ min, max 

Vj/ V-j. max.ref 1^2, Nga , (as in Table 3) 

= 270.0° C 
T””= 220.0° C 

Nga= 10 
Np = 80 

Natr=7 

N,i„,= 100 
q=15 
a = 2 
Po = 0.99 
Pm = 0.001 

Ng,max 20 

Nchr = (10+1) X 7 = 77 
Wi = W2 = 0.25 X 10® 


Table 4; Computational Parameters Used in the Reference Run 


single objective function optimizations (i.e., by solving Eq. 6 using only one of the tvvo 
objectives and dropping the other objective function completely). These two limiting 
solutions show up as the straight lines in Figs. 9-11. The intersection of these lines gives 
the utopia. Point U; is unattainable, but represents an ideal. The preferred solution is often 
taken^-’ to be the point nearest to utopia. These points are denoted by O; in Rgs. 9-11. 
The subscript, i , on points U; and 0; corresponds to the integer in the value of the initial 
water concentration, [W]o . It must be emphasized that the geometrical construction 
required to obtain points Oi will depend on the scales being used in plotting the Paretos. 
However, this is the best one can do in absence of additional information. 

Each point on the Pareto sets in Figs. 9-11 corresponds to a different VxA^T^naxjef 
history as well as Tj. Three of these histories (for tf / tf^f = 0.45, 0.54, 0.65 on the Pareto 
for the reference case of [W]o = 3.45%) are shown in Fig. 4. Fig. 12 shows the optimal 
values of Tj corresponding to the Paretos for all three values of fWjo . 

In Fig. 13, the computed Pareto sets are compared with the Pareto sets obtained by 
Sareen and Gupta^® for one case. The Paretos obtained in the present case indicate superior 
performance compared to that obtained earlier. This indicates that ‘freezing’ of the shape 
of the pressure history during optimization does not give the best optimal solutions, since 
the search space gets limited by this procedure, a drawback that does not exist in the 
current, more general, study. 

Fig. 14, shows the optimal history of VTA^T,iiia.x,icf for the preferred solution, O3, and 
compares it with values currently in use.^^ It is observed that the vapor rdease rate is 
finite at the beginning in the case of the optimal solution, in contrast to the present 
operation. Also, the maximum value of the optimal vapor release rate is lower. Fig. 15 
shows that higher vapor release rates are required for higher values of [W]o - Fig. 2 shows 
the variation of the dimensionless pressure, FI » with dimensionless time for the preferred 



solutions (O2 - O4). Much lower pressures are required as compared to the current 
operation. Fig. 16 shows the variation of the dimensionless temperature with dimensionless 
time for the three preferred solutions. The rise of the temperatures are observed to 
increase more slowly as compared to current values. This is because of nonzero vapor 
release rates at the begirming (leading to higher amounts of vaporization). In Fig. 17, the 
number average chain length, Pn , vx the dimensionless time, is plotted for the three 
preferred solutions. The plateau observed currently is found to be missing in the case of 
the optimal runs. This is similar to the observation made by Sareen and Gupta.^ In Fig. 18, 
the monomer conversion, conv, is plotted against the dimensionless time. Delayed 
attainment of the equilibrium \'alues is observed. In Fig. 19, [C 2 ]/[C 2 ]£r <£5 (where subscript 
3 indicates that the reference value, [C 2 ]£ref, for [W]o = 3.45% is used for all three values of 
[W]o) is plotted against the dimensionless time for the preferred solutions as well as for 
the current operation. Use of the same denominator for all three [W]o helps compare the 
values of [C 2 ] better. It is observed that the amount of cyclics are much smaller than 
produced currently, as also indicated in Fig. 13. The quality of the product produced under 
optimal conditions is, thus, observed to be far better (while the total reaction time, tf, is 
also substantially reduced leading to higher plant capacities). 

The effect of varying the several computational parameters given in Table 3 and 4 
on the Pareto solution for [Wjg = 3.45% was then studied. If the population size, Np , is 
decreased from the reference value of 80 to 60, the same Pareto is obtained. li\e could, as 
well, choose to use N? = 60 since the number of function evaluation would be lower. 
However, increase in Np to 100 led to severe memory and CPU time requirements, and 
only three generations could be studied. A few extra feasible solutions were found for Np — 
100 than shown in Figs. 5 - 8. The evolving Pareto appeared to be identical. On increasing 
the value of Ha- from 7 to 10, we obtained Pareto solutions which differed onK- in the fifth 
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Fig. 14 ; Variation of the dimensionless vapor release rate with dimenaonless time for the 
current and optimal (preferred solution, O 3 ) cases; [W]o = 3.45% . 
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Fig. 15 : Dimensionless vapor release rate histories for the preferred solutions for [W]o = 
2 . 52 %, 3 . 45 %, 4 . 43 % . 
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to seventh decimal places. This increased accuracy is, thus, not necessary and the reference 
value Nstr = 7 is alright. Similarly, increase in the value of Nga from 10 to 16 led to no 
significant changes in the Pareto set. The effect of varying the value of the crossover 
probability, po, from 0.99 to 0.77, again, led to no significant changes in the optimal 
solutions. The effect of a change in the value of the mutation probability, pm, from 0.001 to 
0.1, led to fewer points on the Pareto set (but the ones obtained were identical to those in 
the reference case), but convergence was attained more slowly (at the 18th. generation). 
Also, higher pm led to some erratic behavior in the feasible points obtained in each 
generation in the beginning (this was probably the cause of slow convergence). Decrease 
as well as minor increase of the window provided for VTA^T,maxjef (to factors of 0.2 and 1 .2 
or 0.7 and 1.7, of the values in Ref 28) led to fewer points on the Pareto (through the 
same values were obtained for the points). The reference values of the computational 
parameters in Tables 3 and 4, thus, appear justified. 
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Chapter 4 


Conclusions 

Pareto optimal solutions have been obtained for the multiobjective optimization 
problem described in Eq. 6. The vapor release rate history, ViCt), and the value of the 
jacket fluid temperature ,Tj, have been used as the control (optimizing) variables. An 
adapted NSGA (Nondominated Sorting Genetic Algorithm) technique has been used to 
obtain the solutions. Use of this niche based search technique was found to be superior to 
the use of Pontryagin’s minimum principle for dynamic optimization.^ In fact, our attempt 
to extend our earlier work^ using Pontryagin’s principle (for single objective optimization 
of the semibatch nylon 6 reactor) to the more difficult problem described in Eq. 6, failed 
due to severe numerical problems. This indicates the superiority of the present technique. 
It is found that there is considerable scope for improvement in the operation of the 
industrial reactor (reduced total reaction times, reduced cyclic formation, while 
simultaneously maintaining the values of monomer conversion and the number average 
molecular weight). We understand that some improvements have already taken place on 
the industrial reactor recently, along the lines indicated in the present work. This work has 
been partly funded by a grant fi-om the Research Centre, Gujrat State Fertilizers Co. Ltd., 
Vadodara, India. 
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Chapter 5 

Suggestions For Future Work 

Though this work has predicted considerable improvement in the performance of 
the industrial semibatch nylon 6 reactor, a more fundamental optimization study is called 
for. Industrial implementation of optimum vapor release rate histories or pressure histories, 
predicted by our work may not be an easy task, if the industrial reactor studied here is not 
capable of taking care of very sophisticated control. 

It is well known that good models are a prerequisite for any optimization study. In 
this study, the activity coefficients of the monomer and water used were those obtained'^ 
by curve-fitting plant-data under one set of operating conditions. It is necessary to re- 
confirm whether the coefficients used in the expressions of yn, and yw are applicable under 
the different conditions corresponding to the optimal solutions. Use of more fundamental 
relationships for the activity coefficients could also to be studied. 
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Appendix 1 


Al.l Mass and Energy Balance Equations ^ 


^ = -k,[C,][W] + k,'[s.] - k 3 [c,l 4 o + k.V. -[s,]) - 


dt 

dt 


= k,[c,]lwl - k;[s,] - 2k3[s,]u. + 2k;[Wl(u, -[s,]) - k3[s,][c,] + k/iS;] - k,[sJ[C3] 


+ks[S3]+[Si| 


0.U3R3„+0,018R.„ 


dUp 

dt 


= k,[c,]lwl - k,’[S,] - k3H= + k,jw](n, - n,) + k 3 [Wl[C 3 ] - k/ls^l + n, 


0.113R,„-0.018R 

\TTi .1 


dh 

dt 


k, [c,][W] - k '[s,] + k3[c,])ao - k3'(|^„ - [s,]) + 2k3[C3j!i3 - 2 k'(^, - [s,] - [S,. 


+2kJW][c,]-2k/[s,] + ^t, 


0.113R_+0.018R„ 


dU; 

dt 


= k,[C,][W] - k,’[s,] + 2 k,nJ + ik;[W](u, - n,) + k3[C,](u, + 2 ti,) + k,'(n, - 2 yl, -[s 


0.113R„ +0.018R. 


'^^5[^:](f^0 ^\) + “^^5 (M-O ~ (■’•1 [^ 2 ]) 4kj[W][C3] ~ [^ 2 ] M-2 

^ = -k^lCjllW] + k;[S 3 ] - k,[C 3 ](i, + kjV, - [S,] -[S 3 ]) ^[C,.] J 3 

= -k, [C,][Wl + k,'[s,] + k 3 n; - k;'! Wl(n, - U 3 ) - k 3 [C 3 ][W] + k/IS;] 

, , ^O.lOR^^+O.llSR,^. 

-R,^./F + [W] ^ : 


0.113R^+0.0i8R,^ 


dt 


— = -(0.11 3R,^+ 0.01 8R,^) 


contd. Alb 



dt 


= UA X-T + 


1000 1 


ir.(-AH,) -[o.n3R,„C;„ + 0.018R,.,.,C: 


x(T-T,) +C;,„Jo.1 13R^ +0.018R,^](T-T,)) x { 


,muc 


+2.0925xl0'’(T-T,j 



d[M'] R,„ V.lM’l 


dt Vg Vg([M'l+[W''‘ 

l+M) 


dlW'] _ R„ VJW| 

dt V, V,([M'1 +[W'-]+[n-]) 


d[N^'] VjN^j 

dt " Vg([M'']+[W'']+[N"j) 



vm 



dt 


= Vx 


Closure Conditions: 


[S,! = [S:] = [S,] 


l-ti ) 

1 ^ 1 1^0 


contd. Ale 



A1.2 Correlations Used For Mndellinq ^^^^- 


Expression for Vt 

Stage /; Vt - 0 

Stages 2-5-. Vr=R,^+R^ 


y.fdp^ _ [\r]+[w^]-f[N^]fdT'; 


RTVdt; 


+ V„ 


T 


dt; 


Rates of vaporization 

R., = F(k,„a)^([c.]-lc,i;) [c,1;=1M[m1RT 

R™ = F(k,,„a) j[Wl -[Wi;) +F(k,,a)J[Wl -[Wli 

_ [C,l(p-/„Pr) 

‘ (T„pr-p) 


lwi: = ^-;^[W]RT 

{ w 


Pressure 


P = {[M"]+[W'']+[N"]jRT 


Equations for activity coefficients 

f[c,1 


monomer conversion = 1.0 

P™r-P 


Tm - PmO +■ 

Tw = PwO + 


— — ^°-x(monomer conversion) 


0.95 

^'"'o 95 x(nionomer conversion) 


Vapor pressure 


7024.023 


ln[pr(kPa)/l01.3] = 13.0063-—^ 


3816.44 


ln[pr(kPa)/l01.3] = 


Diffusion coefficients 



Latent heats of vaporization 
Tr = 473.15 K 

= 34.2559 kJmor' 

= 51.0193 kJmor' 
C;„ = 1.6426 kJ kg-' K-' 

= 1.9963 kJ kg-' K'' 


Heat transfer coefficient 


hi (kJ h-' m'^ K-') = 


[ri(Pa_sec) 


0.17 


1 thickness(m) 

hi ks3 

Correlations for mixture physical properties (liquid) 

CU.x (hJ kg-' K-') = 2.0925 + 2.0925 X 10-' [T(K) -273.15] 

p (kg m-') = 1000{l.l238 - 0.5663 x 10-'[T(K) - 273.15]} 

k (kJ h-' m-' K-') = 0.7558 

Correlations for mass transfer coefficients 

6ri( poise) 

360ti( poise) 

^Sc.i 


p®. 


contd. Ale 
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n) Quiescent (q) desorption (nb<P) 



(ki.wa)^=(k,,„a)^_^ 

(k,.„a)^ = 0 


[for i = m or w] 


(ii) Bubbly (b) desorption (Qb>P) 

[w]-K]' 

Kf 

a = 1.81N;",“ 


If (j <0^ : 

(k,,a)Jh-') = 677xlO-‘Nj,V” X3600 
I>,..a),=(k,,„a)^^ 


If <j>o^ : . 

(k, ^a) jh-’) = 2.45 X X 3600 

log,o(l) = 522(G-a,)N^/' 

(ki.wa).. =(k,^,a) (j) 



>6 


contd. A If 



Viscosity correlation (neglecting effect of water) 
ri^(cp) = 2.7969 x lO"" exp[3636.364/T(K)] 


for T> 473.15 K 


H 


100 kg mixture 
, kg polymer 




.. n-nmi 


= lim-^ = lim 

c->.o c c-'O ri„ 


1875.0 
L~ T(K) 


+ 4.678 


M„ 


\0 75 


J k5424; 


M„ = 113(a,/iio 
= llSju^/p.! 

c(kg polymer/100 kg mixture) = 1 1.3|Ui 


For c < 21.0 /[ti] 


~ = 1.0- 0.3 102c[r(l + 0.0575 (c[ti])' - 0.525 x 10'-(c[ril)' + 0.2305 x 10 ''(c[ti1)" 
For c > 21.0/[ri] 

log,oTi(poise) = 51og,o(c,'Mr) - 12-3097 

/or > 315 and Mw> 5000 

else 

log,. t^poise) = log,.(c;M,)- 3.503 
where 

Cj(gm polymer/cm^ mixture) = 11.3 x 1 O'" pu, 


- 0.3663 X 




C HS~)^ jr)~ (yjj-j,. mUL 



